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Abstract. We present a new method to estimate the Hubble constant Ho from the measured time delays in 
quadruply imaged gravitational lens systems. We show how it is possible to get an estimate of Ho without the 
need to completely reconstruct the lensing potential thus avoiding any a priori hypothesis on the expression of the 
galaxy lens model. Our method only needs to assume that the lens potential may be expressed as r°‘F{9), whatever 
the shape function F{9) is, and it is thus able to fully explore the degeneracy in the mass models taking also into 
account the presence of an external shear. We test the method on simulated cases and show that it does work well 
in recovering the correct value of the slope a of the radial profile and of the Hubble constant Hq - Then, we apply the 
same method to the real quadruple lenses PG1115+080 and B1422+231 obtaining Hq = km s“^ Mpc“^ (68 
% GL). 

Key words, gravitational lensing - cosmology : distance scale - quasars: individual: B1422+231 - quasars: indi¬ 
vidual: PG1115-I-080 


1. Introduction 


Most of the ways of measuring the Hubble constant Hq 
involve a form of distance ladder, which utilizes a num¬ 
ber of astrophysical standard candles and standard ruler 
relations, and are calibrated locally by a geometrical tech¬ 
nique such as trigonometric or dynamical parallaxes (e.g., 
Madore et ah, 1998). A few methods involve no distance 
ladder: good examples are (i) inferring the distance of 
the SNell from their light curves and spectra by mod¬ 
elling their expanded photosphere ( ^chmidt et ah, 19^ ) 
and (ii) comparing the Hq - independent angular extent of 
galaxy clusters to their Hq - dependent depth as deduced 
by the X-ray emission and the Sunyaev-Zeldovich mi¬ 
crowave background decrement due to the cluster itself 
( Hughes fc Birkinshaw, 1998 ). 

But the most promising one step method may be con¬ 
sidered the one first proposed by Refsdal in 1964 and be¬ 
came feasible only recently thanks to the now available 
instrumentation. The principle of the Refsdal method is 
quite simple. If a QSO is multiply imaged by the gravi¬ 
tational lensing effect of a galaxy along the line of sight, 
the light rays coming from the different images follow dif¬ 
ferent optical path and thus arrive to the observer with a 
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time delay among each other. It is easy to show that these 
time delays are proportional to the inverse of the Hubble 
constant Hq and to a factor which depends only on the 
lensing potential and the source coordinates. Having mea¬ 
sured the time delays and estimated the lens dependent 
factor by the images configuration, one can then obtain 
a direct estimate of the Hubble constant avoiding all the 
problems and possible systematic errors connected to the 
distance ladders. There are actually more than fifty multi¬ 
ply imaged systems (both double and quadruple) and the 
number of them with measured time delays is increasing 
( Bchechter, 2000| ) so that the prospects of obtaining an 
accurate estimate of Hq from gravitational lenses is quite 
good (Koopmans, 2001). 


However, there is still a major problem connected to 
the modelling of the lensing galaxies since there are often 
different models which predict the same images configu¬ 
ration and the other lensing observables. Thus lens mod¬ 
elling is the major source of uncertainty in the Refsdal 
method. There are two ways to compensate for our lack 
of knowledge about the lens galaxy. The first is assuming 
an exact parametric form for the galaxy model and then 
determine its parameters by fitting to the images posi¬ 
tions and time delays ratios (and, eventually, to the flux 
ratios). However, it is clear that this approach strongly un- 
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derestimates the uncertainty connected to the lens mod¬ 
elling. Even if one is able to find a parametric galaxy 
model which is dynamically possible and which repro¬ 
duces the image properties with acceptably low one 
still has to aggressively explore all other classes of models 
to get the true uncertainty on Hq. But this is not pos¬ 
sible with the parametric approach since there is an un¬ 
avoidable limit to the number of parameters which can 
be used to describe the galaxy fixed by the number of 
available constraints. Thus parametric techniques may ex¬ 
plore only the simplest models, i.e. one is restricted to 
a narrow area in the space of the models. To explore it 
in a systematic fashion one has to follow the second ap¬ 
proach using a representation of the galaxy which is as 
general as possible and thus not restricted to a particu¬ 
lar form. One possibility in this sense is to pixelate the 
galaxy map and consider each pixel as an independent 
mass element. This is what is done in the pixellated lens 
method ( $aha fc Williams, 1997 ; Williams fc Saha, 2000 ) 
which indeed is a valide alternative to the usual parametric 
techniques. Introducing as less as possible constraints on 
the reconstructed mass distributions, the pixellated lens 
method is very efficient in exploring the models space, 
but it has also the risk of overestimating the uncertain¬ 
ties on Hq since one has almost no way to control if the 
reconstructed models are physically meaningfull or not. 
To overcome this difficulty we have elaborated a new ap¬ 
proach which aims at estimating the Hubble constant from 
a detailed exploration of that region of the models space 
which is compatible with the lensing observables and, at 
the same time, is physically well motivated. To do this we 
have to lose some generality since we introduce a lensing 
potential which has a well defined radial profile, but we 
still do no hypotheses on the angular part and take also 
into account the presence of the external shear. Even if 
less general than the pixellated lens method, our semi- 
analytical technique may be considered as a compromise 
between the usual parametric technique and the full non - 
parametric approach. We do not introduce any defined 
galaxy lens model, but we take into account all the models 
which give rise to a lensing potential of the form r°‘F{9), 
whatever the shape function F{9) is, and finally select only 
the ones which reproduces the lensing observables and are 
physically well motivated. Thus our method is still able to 
carefully explore the models space and, at the same time, 
it does not introduce any overestimate connected to the 
inclusion of unphysical models. 

The plan of the paper is as follows. In Sect. 2 we 
do some general considerations on the lens equations and 
write down the relations which will be used to build the 
method. This is presented in Sect. 3 where we will show 
how it is possible to algebrically manage the lens equa¬ 
tions to finally get a system which is numerically solvable. 
Since the set of equations is nonlinear, we have elabo¬ 
rated a simple algorithm that allows us to recover all the 
solutions and exclude the unphysical ones. The code and 
the selection criteria are then presented in Sect. 4, while 
the following Sect. 5 is devoted to the application of the 


method to simulated cases and to the description of how 
we extract the Hubble constant estimate from the set of so¬ 
lutions. Having so checked that the method indeed works, 
we can now apply it to real systems; this is the subject of 
Sect. 6 where we discuss the real lenses PG1115-I-080 and 
B1422-I-231 and obtain our final estimate of Hq. Sect. 7 is 
then devoted to conclusions. 


2. Lens equations and time delays 

Let us begin with some general considerations on the esti¬ 
mate of Hq from the time delays between two different im¬ 
ages of the same source. To this aim let us fix a cartesian 
coordinate system with origin on the lens galaxy centre 
and (cc, y) axes along the main axes of the galaxy itself, 
positevely oriented towards West and North respectively. 
Let (r, 9) be the polar coordinates with the 9 angle mea¬ 
sured counterclockwise from North. Thus we have: 


X = r sm t 


y = r cos 9 . 


The time delay of a light ray deflected by the galaxy lens¬ 
ing effect is given by (Zhao & Pronk, 2001): 


At{r, 9) = h Vioo 


- rxs cos {9 - 9s) + 


+ 


( 1 ) 


being (r,9) the image position, (rs,9s) the unknown 
source position and 9) the lensing potential. In Eq.(P, 
h is the Hubble constant Hq in units of 100 km s“^ Mpc“^ 
(i.e. Hq = lOO/i km s“^ Mpc“^), while tiqo is a typical 
time delay estimated for a given set of cosmological param¬ 
eters (r2m,HA,flfc) assuming Hq = 100 km s“^ Mpc“^. 
The typical time delay is defined as: 


TlOO 


_ / DolDos \ (1 -I- zl) 
V Dls ) c 


( 2 ) 


with the usual meaning for the angular diameter distances 
Dol,Dos, Dls', Zl is the redshift of the lens. 

From Eq.(||) it is easy to get the time delay between 
two images with coordinates {ri,9i), (rj,9j) respectively: 


Atij = Ati — Atj = h ^Tioo X 


-i;{ri,9^)+il}{rj,9j)} . (3) 

According to the Fermat principle, the images lie at the 
minima of At, so that the lens equations may be simply 
obtained minimizing At. We get: 


d 


dr 

i_a 


At = 0 




r - Ts cos {9 -9s) = , 


( 4 ) 


At = 0 


Ts sin {9 -9s) = 


1 d'lp 


( 5 ) 
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The lensing potential Q) niay be splitted in two terms 
and it is usually written as: 


-Ipir, 0) = ijjlens (r, 0) + i! shear (f, 0) ■ 


( 6 ) 


'^^4’lens = 2K(r, 9) 


(7) 


t(r, 9) = 


nr,9) 


4’lens = r°^F{9) . 

Introducing this expression into Eq. ( 0 ) we get: 


nr, 9) = -r° 


n{9) 


(PF{9) 

d9^ 


( 8 ) 


(9) 


As an example, we remember the isothermal ellipsoidal 
model (Blandford & Kochanek, 1987) which is obtanined 
from Eq. (H) by putting : 


The first one is connected with the galaxy lens models 
through the 2D Poisson equation: 


= 1 ;F{9) (X \J cos^ 9 + q~'^ sin^ 9 . 


being K(r, 9) the convergence, i.e. the adimensional surface 
mass density defined as ( Schneider et al., 1992 ): 


These considerations lead us to be confident that our 
assumption on ipiens is quite general and allow us to ex¬ 
plore in detail a wide area in the space of models. 

Let us turn now to the second term in Eq.(^. This is 
the potential of the external shear which is introduced to 
take into account (to the lowest order) of the presence of 
other galaxies of the group which the lens galaxy belongs 
to. We have: 


with T,crit = {nnT^G)(DosID olDls)- It is not un- 
usual to suppose that ipiens may be factorized (see, e.g., 
Witt et al., 2000); motivated by this, we will assume that: 


'^shear — 2^^ COS (2^ 2^^) 


( 10 ) 


Eq.(^ tells us that our hypothesis on ipiens simply means 
that the mass models we are considering are cuspy power - 
law and may be spherical or not according to the expres¬ 
sion of the shape function F{9); since many galactic mod¬ 
els predict a single power - law surface mass density or may 
be approximated (at least in the external part which is the 
one we are mainly interested in) by Eq.(^, we are confi¬ 
dent that our hypothesis on ipiens is not too severe and is 
quite enough general to explore the mass models degener¬ 
acy. 

It is worthwile to point out that lensing potentials of 
the kind in Eq.(^ have yet been taken into consideration 
by different authors in literature. Witt, Mao & Keeton 
( 2000 ) used the same potential as in Eq.(|^ to analyse the 
effect of changing the slope a of the radial profile and the 
shear strength 7 (see later) on the estimate of the Hubble 
constant from the time delays measurements without do¬ 
ing any hypotheses on the shape function F{9). On the 
other hand, introducing a general but parametric expres¬ 
sion for F{9), Zhao & Pronk (2001) were able to develop 
a semianalytical technique to reconstruct the lensing po¬ 
tential in quadruply imaged systems having assigned the 
boxiness parameter. Their work has been then general¬ 
ized by Cardone et al. (2001) which have implemented 
another semianalytical method which is able to recover 
all the lensing potential parameters without assigning ab 
initio any of them. More recently, Evans & Witt (2001; 
see also Hunter & Evans, 2001) have studied the gravita¬ 
tional lensing properties (i.e. number of images and flux 
ratios) of the family of scale-free galaxies with flat rota¬ 
tion curves. As the authors note, the lensing potential of 
these models are of the form in Eq.(^) with a = 1 and F{9) 
given. Finally, we also note that there are other popular 
parametric lens models which may be reduced to Eq.(&. 


being ( 7 , 9^) the strength and the position angle (i.e. the 
angle between the direction of the shear and the main axis 
of the lens galaxy) of the shear itself treated as a pseudo - 
vector. 

Introducing Eqs.(||) and (0) into Eqs.(|l) and d) the 
lens equations become: 

r — Ts cos {9 — 9s) = ar'^~^F{9) — yr cos ( 2 d — 2 d-y) , ( 11 ) 
Ts sin [9 — 9s) = r°‘~^ fi{9)F{9) + yrsin (2d — 29^) (12) 


having defined: 

1 dF{9) 


im = 


F{9) d9 


(13) 


Note that this latter dehnition implicitly assumes that the 
shape function F{9) never vanishes in the positions of the 
images. Thus the method we develop will work only for 
those potentials satisfying this constraint; however, this 
is not a seriuos problem as it may be shown analysing 
many of the potentials in literature. We will also assume 
that /i(d) never vanishes in the positions of the images; 
this constrains us to consider only not spherically sym¬ 
metric potentials since for the spherical potentials F{9) 
is constant and thus /i(d) identicaly vanishes. On the 
other hand, lensing galaxies are usually elliptical ones and 
there are also diffe rent evidences that dark halos are not 
spherical (see, e.g., ^ackett, 1999 and references therein). 
So these considerations suggest us that to consider only 
not spherically simmetric potentials is not a serious limi¬ 
tation. It is useful to solve Eq. ( 0 ) with respect to ipiens', 
we simply get: 


Ipiens iri,9i) = rfF{9i) = 

rl - riTs cos (di - 9s) + 'frj cos {29i - 29^) 


(14) 


where hereinafter we add a label (running from i to 1) to 
the image coordinates to distinguish among them. Note 
that this relation diverges if a = 0 , but this hypothesis 
may be excluded since 0 = 0 means that ipiens does not 
depend on r which is very unlikely. Introducing Eqs. (10) 
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and ( PH ) into Eq. (i and the result into Eq. i) , it is then 
simple to show that the time delay between two images i 
and j is: 

Aty = At{ri,ei) - At{rj,9j) = AU - Atj = 


''"100 


2a 


{(a - 2){rf - Tj) + 2(1 - a)rs x 


[n cos {9i - 9s) - Tj cos {9j - 0s)] + (a - 2) x 

[r^cos(29i — 29j)—rjCOs(29j — 29-y)]} . (15) 

Eq.(p5D shows us that to estimate the Hubble constant h 
from the measured time delay between images i and j one 
has to know the slope a of the radial profile, the source co¬ 
ordinates (rg, 9s) and the shear parameters (7, 9^). It also 
shows us that it is not strictly necessary to completely re¬ 
construct the lensing potential ipiens i-e. we do not need 
to know anything about the shape function F(9) if we 
were able to find the set of quantities (a, Ts, 0s, 7, ^*7)- To 
this aim one usually adopt a parametric approach giv¬ 
ing a lens mass model which depend on some unknown 
parameters. These latter, the source coordinates and the 
shear parameters are then found out using a minimiza¬ 
tion technique which aims at recovering those values of 
the parameters such that the predicted images positions 
and time delays agree well with what is observed. This 
approach is plagued by the contrast between the need to 
introduce many parameters to build up a reliable and ac¬ 
curate galaxy model and the fixed number of constraints 
(eigth from the four image positions and two from the 
time delay ratios). As a consequence the number of de¬ 
grees of freedom may be too low to get an accurate 
enough estimate of the parameters. To increase the num¬ 
ber of degrees of freedom one may also consider the flux 
ratios, but this is dangerous si nce these latter may be 
contaminated by microlensin g ( Chang fc Refsdal, 1979 , 
Koopmans fc de Bruyn, 2000 ) and other secondary ef- 
f ect, such as, e.g., su bstructure in the lens galaxy 
( Mao fc Schneider, 199^ , which are difficult to handle. 
Beside, there are often different models which fit well the 
same set of observations so that the parametric approach 
is unable to fully explore the degeneracy in the lensing 
potential leading to an underestimate of the systematic 
errors on h. In the next Section we will show how it is 
possible to overcome this problem and obtain an estimate 
of the Hubble constant taking into account the degeneracy 
in the lens mass models. 


3. A new semianalytical method to estimate Ho 


Eq.(|l5D tells us that it is possible to estimate the (di¬ 
mensionless) Hubble cosntant h having measured the time 
delays between any pair of images and their positions pro¬ 
vided that one is able to get the slope a of the lens po¬ 
tential radial profile, the source coordinates {rs,9s) and 
the shear parameters (7,0..^). Note that it is not strictly 
needed to know the exact expression of the shape function 


F{9) if we were able to recover the above parameters. This 
suggests that it should be possible to work out a method 
which leads to an estimate of h whatever is the shape func¬ 
tion itself provided that the lensing potential may still be 
written as in Eq.(^ with ipiens as in Eq.(|^). Actually this 
is true as we will show in this Section. 

As a preliminary consideration, let us observe that the 
lens equations may be written for each one of 

the images thus giving us eigth linearly independent equa¬ 
tions. Whatever is the way we manipulate these equations 
we may still obtain the same number of independent equa¬ 
tions which may be used to determine (at least in prin¬ 
ciple) eight unknown variables. In these equations F{9i), 
F{9j), F{9k) and F{9i) are not functions, but numbers 
which we may rename as {Fi,Fj,Fk,Fl)] a similar dis¬ 
cussion also holds for the values of /i(0). So we have to 
determine these set of numbers and not the shape func¬ 
tion itself. The degeneracy in the lensing potential may be 
reformulated as follows : all the shape functions F{9) such 
that their values in the positions of the images are equal 
to the set of numbers {Fi, Fj, Fk, FI) are acceptable since 
the predicted image positions and time delay ratios are the 
same as the observed ones. This observation explains why 
we do not need to give an a priori expression for the po¬ 
tential ipiens '■ we do not solve for the potential itself, but 
only for the values of the potential in the positions of the 
images. These latter are all what we need to finally get an 
estimate of the Hubble constant. As we will see, however, 
we really do not solve with respect to {Fi, Fj, Fk, FI), but 
with respect to (/li, /Ij, flk, fll), being these latter the 
values of fi (0) in the position of the images. This is related 
to the well known circumstance that the image positions 
are determined by the first derivative of the potential and 
not directly by the potential itself. Our previous discussion 
on the lens models degeneracy still holds when considering 

fiiO)- 

Let us now turn to the bulk of the method. As a fisrt 
step, let us write Eq.(|lH) for the image i and solve it with 
respect to fli. It is easy to get what follows: 

I ___ 

fli rsSm{9^ - 9s) -'yrism{9i - 9j) ' 

Solving Eq. © with respect to r“ ^Fi and inserting 
this solution in the previous relation we get the following 
equation: 

a _ _ ri-rs cos (0^ - 9s) + 7^ cos (20^ - 29^) 
fli Ts sin (0i - 0s) - 7ri sin (0i - 0^) 

where we have introduced a new variable 7^. Eq. 
may be rewritten for each one of other three images us¬ 
ing the right coordinates. Thus we get a system of four 
independent equations in eight unknown variables, i.e. 
the four quantities {r]i,rij,r]k,rii), the source coordinates 
(rs, 0s) and the shear parameters (7, 0-y). To close the sys¬ 
tem we need other four independent equations. To this 
aim, let us turn back to Eq. dH) written for image i and 
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solve it with respect to r“ ^Fi. Inserting the solution into 
Eq.® after some simple algebra one finally gets: 


4. Solving the system of equations and estimating 
the Hubble constant Hq 


[1 + 77 ^ tan ( 0 * - 6 »s)]rs cos ( 0 * - 9^) = 

{! + [! + ? 7 i tan (20^ - 29^)\}^ri cos (20* - 20^) . (17) 

We may write down Eq. (0) also for the other three im¬ 
ages thus obtaining a total of four equations of this kind. 
Adding these relations to the four previously found, we 
finally have a system of eight linearly independent equa¬ 
tions in the eight unknowns {rs,9s), {1,9^), {rii,r]j,r]k,r]i) 
which may be solved (at least numerically). Note that 
there is still a degeneracy in the system since we solve 
with respect to {rjijVjjVkTVi) thus we are not able to 
get the slope a of the radial profile which is the quantity 
we are mainly interested in since it appears in the time 
delay equation together with the source coordinates 
(rg, 9s) and the shear parameters ( 7 , 9^). But we have still 
not used all the informations we have at our disposal. Since 
our aim is to determine the Hubble constant from the time 
delays, it is obvious that we have measured these quan¬ 
tities and so we also know the time delay ratios. Using 
Eq.(|l5|) we may write down the time delay ratios as func¬ 
tion of (a,Tg,0s, 7 ,0^). Having yet found out the source 
coordinates and the shear parameters solving the system, 
this latter relation may then be solved with respect to 
a. Note that now we have all what we need to estimate 
h since we can now solve Eq.dl^) with respect to h itself 
being all the other quantities known. However, we may ob¬ 
tain also the values of lensing potential ip^r, 9) in the posi¬ 
tion of the four images. Actually, having found out a and 
knowing the quantities ( 77 ^, rij,r]k, 77 /), we may immediately 
estimate {fli,flj,flk,fll), then obtain (Fi, Fj, Fk, FI) 
from Eq. 0 ) and finally get: 


7/>(ri, 9i) = r^Fi - cos (20* - 0^) 


(18) 


with: 


Fi = 


Ti - Ts cos ( 0 i - 9s) + 'yri cos {29i - 29^) 


(19) 


To recover the whole set of parameters we need to esti¬ 
mate the Hubble constant from the time delay we use the 
procedure we briefly describe here. First, we substitute 
the f]i from Eq.([l^) into Eq.(|^) and do this for each im¬ 
age finally obtaining a system of only four equations in the 
source coordinates (rg, 0 s) and ( 7 , 9^) that we do not write 
here for sake of shortness. Then, supposing to have solved 
this system, we estimate a from the time delay ratios and 
use this estimate and Eqs. 0) to get (/U,/lj,g/lfc,/lZ). 
Finally, we evaluate {Fi, Fj, Fk, FI) from Eq.(p^) and 
then the Hubble constant from the time delay with the 
values of the lensing potential given by Eq.(p^). Note that 
the most problematic step in this procedure is solving the 
system of four equations in (rg, 0 g) and ( 7 , 9^) since this is 
a nonlinear set of equations and its solutions may be found 
only numerically. This has led us to develop an algorithmic 
to search for the solutions of the system. To avoid intro¬ 
ducing any bias in the search, we give to the algorithm 
Af random starting points for (rg, 0 g, 7 , 0^), where A/” is a 
number fixed by the user^. We then obtain J\f solutions 
which are not all physically acceptable. To select among 
these we have imposed a set of selection criteria: the code 
checks the list of Af solutions and finally retains only the 
ones satisfying the whole set of criteria. Schematically the 
selection criteria we impose and the reasons why we use 
them are described in the following. 

1 . 0 < rg < max(r*,rj,rfc,r;): we impose this cut on rg 
since it is not possible that the source is outside the 
ring delineated by the more distant of the images^. 

2. 0 < a < 2: from Eq.(^) it is evident that the surface 
mass density scales as r““^ so that a < 2 is needed 
in order to have K{r, 9) monotonically decreasing with 
the radius; at the same time, it may be shown that the 
projected mass inside r goes as r“ so that a must be 
positive to be physically reasonable. 

3. 0 < 7 < 'jmax ■ the shear intensity 7 is positive by 
definition and is not expected to be too large since it 
is a perturbation; for this reason we impose an upper 
limit jmax which may be fixed quite arbitrarily; we 
decide to use 'fmax = 0.3, as in Witt et al. (2000), but 
another choice could be to impose 7 = jcrit (being 


The method describeed here thus allows us to get an 
estimate of the Hubble constant h without assigning any 
a priori expression for the shape function and it is thus 
able to fully explore the lens model degneracy whatever is 
the angular structure of the lens which is crucial in deter¬ 
mining the image configurations. Note that to this aim the 
method does not need to use the flux ratios which may be 
seriously affected by microlensing and other systematics 
which may lead to significatively wrong results. Beside, 
the method also takes into account the effect of the other 
galaxies of the group which the lens belongs to through 
the introduction of the external shear whose parameters 
are also obtained. 


^ The code is nothing more but a notebook written for 
MATHEMATIC A, which we have named HERQuLeS {Ho 
Estimate Recovered from Quadruple Lens Systems). 

^ Af should be large to explore a wide region in the parame¬ 
ter space, but not too large to save computer time. The right 
choice must be a compromise between these two different cir¬ 
cumstances. A possible strategy could be to fix A/” = 10000 
and then, eventually, run HERQuLeS more than one time if 
necessary. 

® A nice demonstration of this may be obtained using 
the web tool developped by K. Ratnatunga which generates 
the image of a lens system given the lensing potential, the 
source coordinates and the observational characteristics (see 
^ittp: //mds .phys . cmu.edu/ego_cgi .html). 
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Fig. 1. Histograms of the recovered values of the slope a (right) and of the Hubble constant h (left) for simulation 1 
normalized such that the total area under each histogram is one. 


this latter dehned such that for values of 7 > 'jcrit the 
estimated Hq becomes negative) which depends on the 
particular lens system to be considered (O. Wucknitz, 
priv. comm.). 

4. the shear is approximately well directed: by this we 
mean that the position angle 9-y of the shear must be 
directed towards the supposed origin of the shear itself; 
e.g., if the shear were due to the cluster which the lens 
galaxy belongs to, then 9^ should be aligned with the 
cluster mass center since there are no reasons why it 
should point elsewhere. 

5. {Fi, Fj, Fk, FI) are all positive: this constraint turns 
out from the consideration that the leasing potential 
'4’iens ('c, 0) is positive dehned; since r“ is always posi¬ 
tive, so must be the quantities (Fi, Fj, Fk, FI) in order 
to have a leasing potential physically acceptable. 

6 . the set of parameters so found solve the lens equations: 

this simply means that we insert the solutions in the 
lens equations to check if the solutions found 

is a correct one or the result of a wrong convergence 
of the numerical algorithm. 

7- hmin < h < hjnax '■ this Condition is imposed to avoid 
that the recovered estimate of h is not physical; in 
principle, one should hx {hmin,hmax) = ( 0 , 1 ) to not 
introduce any bias in the estimate of h, but it is also 
well known that different method of estimates never 
give values of h less than 0.3 so that we have used 
[J^rain •! hmax) — (0.25,0.95). 

Obviously, one may also change the order of the con¬ 
straints; we have chosen this one as a tentative to minimize 
the CPU time needed, but nothing prevents to add or re¬ 
move some selection criteria. At the end of the selection 
procedure, one has a set of M. solutions each one with 
attached an estimate of the Hubble constant h. Note that 
these solutions could also be different from each other (and 
actually this is what happens) since the system we have 
initially solved is nonlinear and thus may have more than 
one solution. This is another way to see the lens models 
degeneracy. A final estimate of h may be obtained consid¬ 
ering the mean as central value and as acceptable range 
the one which contains the 90% (or the 68 %) of the values 
so found thus finding out what we will call the 90% CL 
(or 68 % CL) range. As we will see in the next section, this 
will give a quite large range, but this is simply a conse¬ 
quence of having taken fully into account the lens model 
degeneracy. However, we will also see how it is possible to 


narrow the range for h combining the data from different 
lens systems. 


5. Application to simulated systems : refining the 
estimate of Hq 

To test if our method indeed works and to see whether it 
is possible to reduce the uncertainty on the estimate of the 
Hubble constant Hq we have constructed three simulated 
quadruple systems. To this aim, we have used the following 
expression for the shape function: 

F( 6 I) = |1-(5cos(26i-26ip)|^ (20) 


where /3 is a boxiness parameter, (5 is a flattening indi¬ 
cator a nd dp is the position ang le of the lensing pot en- 
tial (see jZhao fc Pronk, 2001 and Cardone et al., 2001 for 
further details). Different choices of the three parameters 
(/3, 5,9p) and of the slope a of the radial profile allows 
to generate different simulated systems provided that one 
have also fixed both the source coordinates {rs,9s) and 
the shear parameters ( 7 , 6 *.^). To compute the time de¬ 
lays one has also to choose the cosmological parameters 
{flm, flfc, h) and the redhsifts of lens {zl) and source 
(zs)- We adopt a flat cosmological scenario for all our sim¬ 
ulations fixing: 


(H^, Ha, Hfc, h) = (0.3,0.7,0.0,0.72) 
and: 

{zl,zs) = (0.310,1.722) —> rioo = 33.37 days arcsec“^ . 

The other parameters were fixed as resumed in the follow¬ 
ing Table 1. 


Id 

irs,es) 

( 7 ,^ 7 ) 

{a, (3, 5, Op) 

51 

52 

53 

(0.09,40°) 
(0.09,40°) 
(0.10,15°) 

(0.20, 36°) 
(0.12,0.36) 
(0.15,66°) 

(1.0,1.0,0.07, 45°) 
(1.2,1.0,0.13,45°) 
(1.1,-0.25,-0.15,32°) 


Table 1. Source coordinates, shear and lensing potential 
parameters for the simulated lenses. 


The results of applying the method to these simulated 
systems is well resumed in Fig. 1 where, as an example, 
we plot the histograms of the recovered values of the slope 
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Fig. 2. Final histogram (normalized such that the total area is one) of the recovered values of the Hubble constant h 
obtained combining the results from the three simulations as described in the text. 


a and of the Hubble constant h for our simulation 1. The 
first interesting result is that the method indeed works 
well in recovering the a parameter which is, perhaps, the 
most important. Considering the 90% range, we get: 

a = for simulation 1 ; 

a = 1 . 18 l'lg' 2 i /o’" simulation 2 ; 

a = 1.13lg'^7 for simulation 3 . 

Comparing these results with the values in Table 1, 
one sees that the method recovers the correct value of a 
with a good enough accuracy since the central value of 
the delineated range for each simulation is very close to 
the input one. The uncertainty is not too large (~ 20%) 
which will impact mostly on the estimate of h. Actually, 
from Eq.([l^) it is clear that a and h are anticorrelated so 
that a little error on a leads to a large uncertainty on h. 
This is clearly shown in the histograms for h which lead 
to the following estimates (90% CL): 

h = 0.76/g'2g for simulation 1 ; 

h — 0 . 74/032 for simulation 2 ; 

h = 0 . 73/030 Z®’" simulation 3 . 


The central value of the range is close to the input h, 
but the uncertainty is quite large (^ 40%). However, if not 
encouraging, this result is not unexpected because of the 
yet discussed anticorrelation between a and h. It is once 
again a consequence of having taken fully into account the 
lens model degeneracy and so it is not a surprise that our 
estimates have a so large uncertainty. On the other hand, 
however, we may try to reduce this uncertainty combin¬ 
ing the results from different lens systems since it is as 
we were marginalizing on the model parameters to retain 
only the Hubble constant. T o combine different resul ts we 
proceed as follows (see also Williams fc Saha, 2000|) . For 


each simulation, we divide the set of values recovered for 
h in bin of 0.01 and to each bin we assign a probability 
defined as the number of values in that bin divided by the 
total number of results. Then we build a combined his¬ 
togram giving to each bin a probability which is equal to 
the product of the three probabilities from each simula¬ 
tion. Finally, we delete from the sample those bins which 
have a combined probability less than 1 % since these val¬ 
ues turn out to be very unlikely. Note that this procedure 
reduces the number of values in the final sample, but we 
believe that this does not lead to exclude any possibly im¬ 
portant value since the excluded bins are really unlikely. 
The result is the histogram in Fig. 2, where we plot bins of 
width 0.1 for sake of clarity. The combined histogram leads 
us to the following estimates for the Hubble constant: 

Hq = 78± 13 km s“^ Mpc“^ at the 68 % CL ; 

Hq = 78/^g km s“^ Mpc“^ at the 90 % CL . 

The uncertainty has been indeed reduced to ~ 20% 
which is a nice result and suggests that it is possible to 
further reduce the range adding other systems. 

It is worthwile to spend some words about the shape 
of the histograms in Figs. 1 and 2. These could suggest an 
estimate of h near ^0.9 where the probability gets higher, 
whilst we report h ^ 0.7-7 0.8 as best estimate. This is due 
to having chosen to use the arithmetic mean as best esti¬ 
mate; even if the last bin has the higher probability, this 
latter is less than the sum of the remaining ones and this 
cause the mean to be lower than what Figs. 1 and 2 could 
suggest. If we had chosen to consider as best estimate of h 
the value corresponding to the bin with the highest proba¬ 
bility, we should get an higher value for h in disagreement 
with the input one. This could be explained if one con¬ 
siders that it is the full ensemble of models which must 
be take into account when determining the Hubble con¬ 
stant and not only the most likely ones. To use the mean 
of the sample as estimate of /i is a simple way to fully 
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Fig. 3. Histograms of the recovered values of the slope a (right) and of the Hubble constant h (left) for PG1115+080 
normalized such that the total area under each histogram is one. 


taken into account the lensing model degeneracy. Our re¬ 
sults from the simulations conhrm this expectation and 
further enforce our adopted procedure for the estimate of 
the Hubble constant. 

We do not discuss here the results obtained for the 
shear parameters and the source coordinates since these 
quantities are less interesting. We limit ourselves to say 
that the central values recovered by the method are always 
close to the input ones, but the uncertainties turn out to 
be quite large. 


6. Application to real systems : PG1115+080 and 
B1422+231 


The simulations described in the previous section have 
shown that the method indeed works recovering the cor¬ 
rect values both of the slope a of the radial profile and of 
the Hubble constant TJq. Given these encouraging results, 
now we apply the method to real quadruply imaged sys¬ 
tems. To this aim, we need a system with four images and 
a good astrometry not only of the images position, but also 
of the lens galaxy centre since this latter will be used as 
the origin of our coordinate system. Beside, we also need 
that there is only one galaxy acting as lens otherwise it is 
not possible to factorize the lensing potential. Finally, to 
estimate the value of the Hubble constant we need that the 
time delays (or, at least, one time delay and the ratio be¬ 
tween any pair of delays) have been measured. Searching 
the CASTLES database ( |Kochanek et al., 2000 ), we have 
found that there are only two systems which satisfy these 
requirement^: PG1115-I-080 and B1422-I-231. In what fol¬ 
lows, we first apply the method to the two systems sep¬ 
arately and then combine the results to get a better es¬ 
timate of the Hubble constant. In all the applications we 
adopt a flat cosmology with (Om,HA) = (0.3, 0.7); the ef¬ 
fect of changing these cosmological parameters is of the 
order of few percent so that we may neglect it. 


There is also a third lens system which may seem to be 
interesting, i.e. B1608-I-656. Actually, this system may not be 
considered since there are two lensing galaxies (probably un¬ 
dergoing a merging event) so that the term ipiens in the lensing 
potential is not separable and thns onr method cannot be ap¬ 
plied. 


6.1. Application to PG1115+080 

The first system we consider is the well known 
PG1115-I-080, first discovered by Weymann et al. (1980) 
and then studied in detail by several authors (see, 


e.g. 

Keeton & Kochanek, 1997 


Williams & Saha, 2000; 

[Zhao 

& Pronk, 2001 

; Cardone et al., 2001 

) with different 


techniques. This system consists of four images (named 
Al, A2, B and C) of a radio quiet QSO at zs = 1.722, 
while the lens is an elliptical galaxy belonging to a group 
of ~ 10 galaxies at zl = 0.310. The center of the group 
is at (vg, 9g) = (20" ± 0.2", —117° ± 3°) and its effect has 
been taken into account as an external shear in previous 
models. We have applied the method to this system us¬ 
ing as images coordinates the ones measured by Impey 
et al. (1998) with HST observations; the time delay be¬ 
tween images Al and A2 is too small to be measured, while 
the one between images B and C is Adsc = 25.0 ± 1.7 d 


(Barkana, 1997) with image B arriving last. There are also 


two different measures of the time delay ratio among im¬ 
ages BC and AB: Schechter et al. (1997) first reported 
fABC — lAtAB/lAtsc = 0.7 ± 0.3, while a later analysis 
by Barkana (1997) found tabc = 1-13 ± 0.18. It is this 
last value which we consider more reliable and use in the 
application of the method. Fig. 3 shows the histograms of 
the recovered values of the slope a of the radial profile and 
of the dimensionless Hubble constant h. Our main results 
(given as 90% GL) are the following ones: 


a = 1.03t°;2o 


h = 0.68t°J5 


( 21 ) 


The slope a of the radial profile is practically equal to 
the one of the isothermal model, consistent with what is 
expected if galaxies are embedded in a dominant isother¬ 
mal dark halo. The central value for h is also very near 
to the previous estimate of the Hubble constant obtained 
for this system. PG115-I-080 has been studied in detail 
by different authors using different techniques. Keeton 
& Kochanek (1997) have used the usual least para¬ 
metric approach modelling the lensing potential as the 
sum of a term due to the galaxy and an external shear 
from the group. They have used various elliptical mod¬ 
els and have found that none of them may fit the im¬ 
age positions without any external shear. A carefull anal¬ 
ysis of the possible systematic effects and of the differ¬ 
ent models also including the external shear from the 
group hnally lead them to estimate the Hubble constant 
as Hq = 5llj^3 km s“^ Mpc“^. Our estimated 90% range 
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Fig. 4. Histograms of the recovered values of the slope a (right) and of the Hubble constant h (left) for B1422+231 
normalized such that the total area under each histogram is one. 


does overlap well with their one which is a very encour¬ 
aging result since our method is completely different. 
We also note that they also considered the isothermal 
model which turns out to fit well the lensing observables. 
Unfortunately their isothermal model is a spherical one 
and thus its lensing potential does not fall in the class of 
the models considered by our method. However, the same 
estimate of Hq is obtained even if one sligthly flattens 
their model. The result is the isothermal ellipsoidal model 
(Blandford & Kochanek, 1987) which belongs to the class 
of potentials (^, as we have yet noted in Sect. 2, and also 
fits well the image positions. Zhao & Pronk have applied 
their semianalitycal method to PG1115-I-080 using a lens¬ 
ing potential which is obtained from Eq.(|^) fixing F{6) as 
in Eq.J^). They examine two different models with fixed 
valued of the boxiness parameter /3 and finally turns out 
with two quite different estimates of Hq ranging from 20 
to 50km s“^ Mpc“^ or from 50 to 90km Mpc“^. The 
first case is marginally consistent with our result, whilst 
the second one is in good agreement. A direct comparison 
with the value of the slope a is not possible since they 
do not report the estimate for this latter parameter, but 
only three vales for what they call the effective power - 
law slope. In an our previous paper (Cardone et ah, 2001) 
we have presented a semianalytical technique to recon¬ 
struct a lensing potential as in Eq.(||) with F{9) as in 
( p0| ) using as constraints only the images positions. Then 
we have applied this method to PG1115-I-080 obtaining 
Fto = km s“^ Mpc“^ still in agreement with the 

present result in Eq. It is also enncouraging that in 
that paper we get a = 1.12 still in agreement with the 
estimate (ill) obtained with the new approach used here. 
Another interesting comparison may be done with the re¬ 
sults ob tained on PG1115-I-080 using the pixellated lens 
method (Williams & Saha, 2000) since this is a not para¬ 
metric technique. Their 90% GL estimate ranges from 30 
to 75 km s“^ Mpc“^ which well overlaps our estimated 
range. Note also that the uncertainties on h from the 
pixellated method are of the same order as our own re¬ 
flecting once again the fact that both methods fully take 
into account the lens models degeneracy. These success- 
full comparisons confirm lead us to be quite confident in 
the validity of our new method since it turns out to give 
results in good agreement with all the previous ones in lit¬ 
erature obtained by using different technique but lensing 
potentials belonging to the same class we consider here. 


Finally we also note that the uncertainty on our estimates 
of a and the wide 90% GL for h are of the same order 
as the ones we have obtained from the simulations thus 
leading us to consider another lens system to reduce the 
estimated range for the Hubble constant. 


6.2. Application to B1422+231 


The second system we consider is B1422-I-231 fisrt discov¬ 
ered by Patnaik et al. (1992) and then observed in detail 
both in optical (Lawrence et ah, 1992; Impey et al., 1996) 
and in radio ( Patnaik et al., 19991 ). The images configu¬ 
ration is quite different from the typical cross which is 
observed in quadruply systems (e.g. PG1115-I-080 and 
G2237-I-030) since we have three images on a side of 
the main lens galaxy and a fourth one very far on the 
other side. This suggest that the source is located near 
a cusp and is not aligned with the galaxy centre as of¬ 
ten in quadruply systems. The four images have redshift 
zs = 3.62 while the main lensing galaxy is at zl = 0.49. 
There are other two galaxies not too far from the main one 
at the same redshift; we will consider them as the sources 
of the external shear and impose that our reconstructed 
shear points in the direction between these two galaxies. 
As input we use the images and galaxies positions de¬ 
termined by Impey et al. (1996) with HST observations, 
while we take as time delays the recently determined val¬ 
ues by Patnaik & Narasimha (2001) even if these latter 
have a very large uncertainty. Fig. 4 shows the histograms 
for a and h obtained applying our method. The main re¬ 
sults we get are: 


a = 1.08t°i3 


h = 0.47t°;3i 


( 22 ) 


The slope a of the radial profile is almost the same 
as for PG1115-I-080 but is determined with a higher un¬ 
certainty. The result on h is qualitatevely different; even 
if the two 90% range do overlap enough, the histogram 
of h values is peaked towards lower ones in clear contrast 
to what has been obtained for PG1115-I-080 and is ex¬ 
pected. Kormann, Schneider & Bartelmann (1994) have 
shown that the lensing observables for B1422-I-231 may 
be reproduced by an isothermal model (a = 1) and an ex¬ 
ternal shear and our result on a is in good agreement with 
their conclusion. This is encouraging since the lensing po¬ 
tential they use is the sume of the isothermal ellipsoidal 
model and and an external shear and thus belongs to the 
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Fig. 5. Final histogram (normalized such that the total area is one) of the recovered values of the Hubble constant h 
obtained combining the results from PG1115+080 and B1422+231 as described in the text. 


class of potentials described by Eqs. (|), d) and dl^) con¬ 
sidered here. This could suggest that the model is not so 
bad, but this conclusion can not be enforced by a compar¬ 
ison with other estimates of Hq simply because this is the 
first time that B1422-I-231 has been used to estimate the 
Hubble constant from time delays. However, we aware the 
reader that our estimate relies on the time delays mea¬ 
sured by Patnaik & Narasimha (2001) which are affected 
by quite high uncertainty. Choosing values for the time de¬ 
lays (consistent within the errors) different from the one 
we have used lead to a different time delay ratio and con- 
sequentely to different estimates of a and h. Given this 
situation, we have decided to still retain the results on h 
coming from B1422-I-231 and combine them with the ones 
from PG115-f080. 

6.3. Combining the results : final estimate of the 
Hubble constant 

Here we combine the histograms obtained for 
PG1115-I-080 and B1422-I-231 to get the final esti¬ 
mate of the Hubble constant. To this aim, we apply the 
same procedure used for the simulations and described in 
Sect. 5. The result we get is shown in Fig. 5 where we 
have used a binning of 0.1 for sake of clarity. Our final 
estimate for the Hubble constant is thus: 

Ho = 581)5 km s'^ Mpc'^ {68%CL) , (23) 

or more conservatively : 

Ho = 58l^^ km s'^ Mpc'^ (90%CL) . (24) 

Note that the error increases significantly when we con¬ 
sider the 90% instead than the 68% CL and this is essen¬ 
tially due to the presence of the peaks to lower values 
due to the unexpected results from B1422-I-231. Should 
this peak be eliminated (considering different time delays 
estimates for B1422-I-231), the 90% range should be re¬ 
duced too. However, also limiting ourserlves to the 68% 


CL range, we note that our final estimate is in good 
agreement with the previous ones in literature. Williams 
& Saha (2000) used the pixellated lens method to es¬ 
timate Ho from quadruple lenses combining the results 
from PG1115-I-080 and B1608-I-656 to finally get Ho = 
61 ± 11 km s“^ Mpc“^ (68% CL) which is in good agree¬ 
ment with our finding. The parametric technique has 
been applied to different lens systems, both double and 
quadruple. We have yet quoted the result obtained by 
Keeton & Kochanek (1997) for PG1115-I-080; here we 
stress that our final estimate of Hq is still in agree¬ 
ment with the value found by them. It is also interest¬ 
ing to compare the result we get with the one obtained 
by Koopmans & Fassnacht (1999); applying a parametric 
technique to B1608-I-656 and considering different mod¬ 
els, they finally quote Hq = 65^6 km s“^ Mpc“^ in satis¬ 
factory agreement with our 68% CL on Hq. Finally, we 
also stress that our result is consistent also with com¬ 
pletely different methods of estimate of the Hubble con¬ 
stant, such as the res ult of the HST Key Pr oject, Hq = 
72±8 km s“^ Mpc“^ ( Freedman et ah, 2001 ), using local 
estimators (as Cepheids and SNIa) and the one from an 
orientation unbiased sample of Sunyaev-Zeldovich clus¬ 
ters, Hq = 651® km s“^ Mpc“^ ( |Jones et ah, 2001 ). 


7. Conclusions 

In this paper we have presented a new semianalytical 
method to estimate the Hubble constant from the mea¬ 
sured time delays in quadruply imaged gravitational lens 
system. Assuming that the galaxy lens potential may be 
splitted into a radial part described by a simple power- 
law profile and an angular part described by a quite gen¬ 
eral shape function, we have been able to write down a 
nonlinear set of equations taking into account also the 
contribution of the external shear. This system may be 
solved numerically allowing us to fully take into account 
the lens models degeneracy; a set of physical constraints 
is then used to select the only reasonable solutions thus 
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avoiding the risk of including in our considerations also 
non physically motivated models. The final set of solu¬ 
tions may be seen as a parametrization of that class of 
lens models which are able to generate the observed im¬ 
ages configuration and the measured time delay ratio. The 
class of models so delineated may be translated in a sam¬ 
ple of values for the Hubble constant Hq which leads us 
to the final estimate of this cosmological parameter. 

To test the method we have applied it to three differ¬ 
ent simulated lens systems varying both the angular and 
the radial part of the lens potential and also the exter¬ 
nal shear parameters. This tests have shown us that the 
method indeed works and it is also very efficient in deter¬ 
mining the slope a of the radial profile with a reasonable 
accuracy. Even if this were not our final aim, the ability 
of our method to find out the a parameter is a very in¬ 
teresting byproduct since this quantity is very usefull in 
modelling the dark halos which may be considered as the 
most important galaxy component responsable of the ob¬ 
served quadruply imaged lens system. The tests have also 
shown us that the method is not efficient in recovering 
the Hubble constant; even if the central value of the 90 % 
range individuated for Hq is very near to the input value, 
the range itself is quite large. This is not an unexpected 
result since we are fully taking into account the lens mod¬ 
els degeneracy and it is well known that this increases the 
uncertainty on the Hubble constant. Motivated by this re¬ 
sult, we have tried to reduce the uncertainty on Hq com¬ 
bining the estimates from different systems. To this aim 
we simply build a combined histograms of the Hq values 
(binned by 0.1 km s“^ Mpc“^) multiplying the probabil¬ 
ities from each system and then excluding from the final 
sample those values of Hq which turn out to have a fi¬ 
nal probability less than 1%. This procedure allows us 
to reduce the uncertainty on Hq by combining the three 
simulated systems. The final estimate of Hq is perfectly 
consistent with the input value used for the simulations 
which confirms us that the method indeed works. 

Given these encouraging results, we have then ap¬ 
plied our method to the real lenses PG1115-I-080 and 
B1422-I-231, which are two of the only three quadruple 
lenses for which the time delay has been measured. As 
regard PG1115-I-080, we find a slope a = l-OStoJo (90% 
GL) which indicates a near isothermal model consistent 
with previous models in literature. The 90% range for Hq 
turns out to be quite large (25 -F 78 km s“^ Mpc“^), but 
we note that the central value {Hq = 68 km s“^ Mpc“^) 
is consistent with the values quoted in literature and ob¬ 
tained with different method. As regard B1422-I-231, we 
find a = 1.08lg'3g ^^0 Hq = 25 -F 78km s“^ Mpc“^, with 
a distribution of values for the Hubble constant peaked 
towards lower ones. As we know, this is the first time that 
this system has been taken in consideration to determine 
the Hubble constant since the time delays have been mea¬ 
sured only recently. Gombining the resulting histograms 


from PG1115-I-080 and B1422-I-231 leads us to the follow¬ 
ing final estimate for the Hubble constant: 

Ho = 58ti^ km s'^ Mpc-^ (68%CL) , 

or more conservatively : 

Ho = 56^26 km s"^ Mpc"^ (90%CL) . 

We note that the uncertainty singnificantly increases 
when passing from the 68% to the 90% range for the esti¬ 
mate of Hq and that the magnitude of the increasement is 
higher than expected when compared with the result from 
our simulations. Analyzing in detail the data, however, it 
is easy to see that this strange behaviour is completely 
due to the histogram obtained for B1422-I-231, which pre¬ 
dicts too many models with low values of Hg. A possible 
explanation of this strange behaviour may be connected 
to the very high uncertainty on the measured time delays 
which translates to an high uncertainty on the recovered 
parameters a and h. Should we have chosen different val¬ 
ues for the time delays (still within the uncertainties), we 
should have obtained lower values for a and higher values 
for Hg thus narrowing the 90% range. Some tests have 
suggested us that this could be actually a possible expla¬ 
nation. However, also considering only the 68% CL, our 
final estimate of Hg is consistent with all the previous es¬ 
timates in literature whatever is the method used and the 
estimators considered. 

The new semianalytical method presented here turns 
out to be a valid complementary alternative to the usual 
parametric approach and to the fully non - parametric 
techniques. Parametric methods are usefull in building up 
galaxy models which may be easily compared to other 
galaxy observables (if possible) also not directly connected 
with the lensing charachteristics. However, they have only 
a modest power in exploring the parameter space since 
one cannot include too many parameters in the model to 
avoid having a number of degrees of freedom too low in the 

minimization. Given that the number of constraints is 
fixed (eight from images positions and two from time de¬ 
lays ratios if the system is a quadruple one), there in an 
unavoidable limit to the accuracy of the galaxy model and 
to the possibility to explore the wide range of lens models 
that may fit the same lensing observables. On the other 
hand, non parametric methods try to introduce as less 
a priori hypotheses as possible thus fully exploring the 
space of the models. However, a fully non parametric ap¬ 
proach as the one adopted in the pixellated lens method 
of Williams & Saha (2000) does not allow to control if the 
models considered are physically motivated or not thus 
risking to overestimate the uncertainties connected to the 
modelling problem. Our new method is less general than 
the pixellated one, but it has the advantage that it se¬ 
lects only models which are physically reliable. Besides, 
combining the results from different quadruple lens sys¬ 
tems helps in reducing the other unidentified sources of 
systematic errors connected to single systems leading to 
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a final estimate for Hq which correctly takes into account 
all the possible sources of errors. 

Further improvements are however still possible. 
Numerical simulations of galaxy formation in different 
cosmological backgrounds have suggested that the dark 
halo d anaity profile may be dcacribed by a simple univer — 
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sallawl (acc, c.g., Navarro et al., 10 9? r Moorc et al., IOOS t- 
Klypin et al., 2001| ). These models lead to leasing poten¬ 
tials which may not be described by Eq.(||) since the radial 
profile is not a single power-law, having different radial 
slopes for the inner and the outer parts. This could suggest 
that the potential we have used is not realistic and thus 
our results on the Hubble constant are wrong. However, 
we have shown that our estimates turn out to be in good 
agreement with all the previous ones obtained both with 
leasing based method and other distance ladders (such as 
Cepheid or Snia). It is worth to note that the leasing ob¬ 
servables are mainly determined by the mass distribution 
in the outer parts of the dark halos and this latter may 
be well described by models with a single power - law ra¬ 
dial profile. Our models differ from the one predicted by 
numerical simulations only in the inner parts which are 
less important in determining the images positions and 
the time delays. It is thus expected that the difference 
does not introduce any serious systematic error. Anyway 
this does not mean that the results from numerical simu¬ 
lations are wrong since the statistic on the leasing systems 
considered is too low. To understand whether these mod¬ 
els are really able to reproduce the leasing observables 
(number and position of images and time delay ratios) in 
quadruply imaged systems and whether they introduce a 
significative change in the estimate of i7o, one has to wait 
for more quadruple lenses with measured time delays. At 
the same time, it should be interesting to generalize our 
method in order to allow also a varying slope a of the ra¬ 
dial profile thus possibly leading also to some constraints 
from the observed quadruply imaged QSOs. 

Finally we note that the number of observed quadruple 
systems to which our method may be applied is consid¬ 
erable so that we have just to wait for the measure of 
the time delays to finally get an accurate estimate of the 
Hubble constant competitive with the ones coming from 
local estimators. 
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